In [1]:
%pylab inline
import pandas as pd
from scipy import stats
import statsmodels.api as sm
import matplotlib.pyplot as plt
from os import listdir
from os.path import isfile, join
import pickle
import itertools
import statsmodels.formula.api as smf
import statsmodels.tsa.stattools
import warnings
from itertools import product
warnings.filterwarnings('ignore')
In [2]:
import statsmodels
statsmodels.__version__
Out[2]:
In [3]:
def invboxcox(y,lmbda):
if lmbda == 0:
return(np.exp(y))
else:
return(np.exp(np.log(lmbda*y+1)/lmbda))
since the task of the third week does not contain any details - lets take half year time period - from january till june 2015 and region №1231 which contains Empire State building. This region is popular popular place.¶
In [2]:
mypath = "H:/Yandex machine learning/finall course coursera/csv_yellow_cab_2015/reg_bin_stat/"
onlyfiles = [f for f in listdir(mypath) if isfile(join(mypath, f))]
print(onlyfiles)
In [3]:
ESB_reg = 1231
month = []
In [4]:
for file in onlyfiles[1:] + [onlyfiles[0]]:
with open(mypath + file, "rb") as f:
region_bined_stat = pickle.load(f)
month.append(region_bined_stat[ESB_reg-1, :])
for i in month:
print(len(i), end=" ")
In [5]:
month_join = list(itertools.chain.from_iterable(month))
time = pd.date_range('2015 Jan 1 00:00:00', periods = len(month_join), freq = 'h')
taxi_data = pd.DataFrame({"taxi_call_num" : month_join}, index = [time])
taxi_data.head()
Out[5]:
In [9]:
print(taxi_data.index.values.shape)
print(taxi_data.taxi_call_num.values.shape)
In [16]:
fig, axes = plt.subplots(figsize = (80, 8))
axes.title.set_size(40)
taxi_data.plot(color="green", title="timeserie", fontsize=25, ax = axes)
frequency = 24
plt.xticks(time[::24])
plt.grid()
#double click on image makes it larger
plot first three month separately¶
In [21]:
#num of hours in each month
month_len = np.array([0, 744, 672, 744, 720, 744, 720])
for i, c in enumerate(["January", "February", "March"]):
taxi_data[int(sum(month_len[:i+1])) : int(sum(month_len[:i+2]))].plot(figsize=(20, 5), title=c, color="green", fontsize=16, grid=True)
26-27 january 2015 - ice storm https://www.nytimes.com/2015/01/27/nyregion/new-york-blizzard.html 8 march - narrow dip may be connected with demonstration http://www.chicagotribune.com/news/nationworld/chi-new-york-city-marks-international-womens-day-with-march-20150308-story.html , or may be just a artifact¶
Dickey-Fuller stationarity test + time-serie STL-decomposition:¶
In [27]:
plt.figure(figsize(60,10))
sm.tsa.seasonal_decompose(taxi_data.taxi_call_num).plot()
print("Dickey-Fuller stationarity test: p=%f" % sm.tsa.stattools.adfuller(taxi_data.taxi_call_num)[1])
#double click on image makes it larger
half year time serie seems to be non-dispersive but it has distinctly visible monthly seasonality. Рypothesis of nonstationarity is rejected with the Dickey-Fuller criterion.¶
plot the same for january month¶
In [28]:
plt.figure(figsize(15,8))
sm.tsa.seasonal_decompose(taxi_data.taxi_call_num[0:744]).plot()
print("Dickey-Fuller stationarity test: p=%f" % sm.tsa.stattools.adfuller(taxi_data.taxi_call_num[0:744])[1])
weekends may be determined from the trend plot and from seasonal graph it is obvious that each day has two maximum values - morning and evening and two minimums - day and night. Night minimum is much deeper.¶
In [10]:
plt.figure(figsize=(15,4))
plt.plot(taxi_data.index, taxi_data.taxi_call_num)
""" moving average """
plt.plot(taxi_data.index, taxi_data.taxi_call_num.rolling(window=24).mean(), 'r')
Out[10]:
From 6-month time serie data it may be concluded that taxi trips have weekly and daily periodicity and also periodicity within one day - night, morning, day, evening. Sunday is most "silent" day then amount of trips starts to grow getting maximum values in the middle of week with a subsequent decrease.¶
Weekly and daily periodicity will be accounted by making a lineat regression model on fourier serie features with periods 7*24 hours (week) и 24 hours (day); morning/evening periodicity should be taken into account by the ARIMA model.¶
T=168 t = 24 W(i) = T/2, T/3, T/4 ... w(i) = t/2, t/3, t/4 ...¶
lets take first 100 garmonics¶
In [13]:
T = 168
t = 28
time = np.arange(0, taxi_data.shape[0])
regression_taxi = taxi_data.copy()
#regression_taxi.loc[regression_taxi["taxi_call_num"] == 0] = 1
#regression_taxi['log_taxi_call_num'], lmbda = stats.boxcox(regression_taxi.taxi_call_num)
In [20]:
#t=28
#tt=48
#ttt = 14
for w in range(1, 100):
regression_taxi["sin_sin_w_" + str(w)] = np.sin(2*np.pi*w*time/T) * np.sin(2*np.pi*w*time/t)
regression_taxi["cos_cos_w_" + str(w)] = np.cos(2*np.pi*w*time/T) * np.cos(2*np.pi*w*time/t)
regression_taxi["sin_cos_w_" + str(w)] = np.sin(2*np.pi*w*time/T) * np.cos(2*np.pi*w*time/t)
regression_taxi["cos_sin_w_" + str(w)] = np.cos(2*np.pi*w*time/T) * np.sin(2*np.pi*w*time/t)
In [21]:
regression_taxi.head()
Out[21]:
In [22]:
regression_taxi.shape
Out[22]:
In [23]:
models = []
for w in range(1, 100):
m1 = smf.ols('taxi_call_num ~ ' + ' + '.join(list(regression_taxi.columns)[1:4*w+1]), data=regression_taxi)
model = m1.fit(cov_type='HC1')
models.append(model)
# fitted.summary(); fitted.params; models[-1].rsquared; params.Intercept
plot r-squared (proportion of the variance in the dependent variable that is predictable from the independent variable) dependence on number of fourier components and similar dependence but for mean squared error of the residuals (sum of squared residuals divided by the residual degrees of freedom)¶
Here we used slighly different features - smth. like "modulation" and take one of the periods not equal to daily period but larger - 28 hours. Obtained result is quite good - up to 90 % of predictable variance.¶
In [27]:
import warnings
warnings.filterwarnings('ignore')
fig1 = figure(figsize=(15,5))
ax1 = fig1.add_subplot(111)
line1 = ax1.plot(range(1, 100), list(map(lambda x: x.rsquared, models)), color="blue", label="r-square")
ylabel("r-squared")
ax2 = fig1.add_subplot(111, sharex=ax1, frameon=False)
line2 = ax2.plot(range(1, 100), list(map(lambda x: x.mse_resid, models)), color="green", label="resid")
ax2.yaxis.tick_right()
ax2.yaxis.set_label_position("right")
ylabel("mse_resid")
legend((line1, line2), ("1", "2"))
plt.xlabel("number of fourier components")
legend()
grid()
show()
plot dependence of correlation of first lag with 168, 24 and 14 lag on number of fourier components¶
In [28]:
statsmodels.tsa.stattools.acf(x=models[3].resid, nlags=168)[168]
Out[28]:
In [30]:
plt.figure(figsize=(15,6))
plt.plot(range(1, 100), [statsmodels.tsa.stattools.acf(x=model.resid, nlags=168)[168] for model in models], color="blue", label='week')
plt.plot(range(1, 100), [statsmodels.tsa.stattools.acf(x=model.resid, nlags=168)[24] for model in models], color="green", label='day')
plt.plot(range(1, 100), [statsmodels.tsa.stattools.acf(x=model.resid, nlags=168)[14] for model in models], color="red", label='14 hours')
plt.legend()
plt.grid()
plt.xlabel("number of fourier components")
plt.ylabel("autocorellation with lag")
Out[30]:
80 components is enought for achiving good performance and quite high decorelation. Фs option is possible to choose 67 components as tradeoff between day and week decorelation¶
In [31]:
models[80].summary()
Out[31]:
function to get the regression prediction¶
In [32]:
def prediction(model, components, feature_num, dataframe):
coefs = np.array([model.params[coef] for coef in model.params.index])
fourier_components = np.hstack((np.array([1]*dataframe.shape[0]).reshape(dataframe.shape[0],1),
dataframe.iloc[:, 1:feature_num*components+1].values))
return coefs.dot(fourier_components.T)
In [33]:
#model = models[40]
pred = prediction(models[80], components=81, feature_num=4, dataframe=regression_taxi)
In [11]:
def prediction_real(df_1, df_2, segment, title, ylabel, size, labels):
plt.figure(figsize=size)
plt.plot(df_1[segment[0]: segment[1]], alpha = 0.7, label = labels[0])
plt.plot(df_2[segment[0]: segment[1]], alpha=0.5, label = labels[1], color="red")
plt.xlabel('time (hours)')
plt.ylabel(ylabel)
plt.title(title)
plt.legend()
plt.grid()
In [56]:
prediction_real(pred, regression_taxi.taxi_call_num.values, (0, 745), "January", 'taxi calls', (15, 5), ("predict", "real"))
In [57]:
prediction_real(pred, regression_taxi.taxi_call_num.values, (0, 4344), "six months", 'taxi calls', (50, 5), ("predict", "real"))
#double click makes plot langer
resids of the model. we make a plot to compare residuals of the model and real taxi values to find is it any significant structure in residual plot.¶
In [60]:
prediction_real(regression_taxi.taxi_call_num.values, models[80].resid.values, (0, 745), "first january week", 'residuals',
(30, 5), ("real values", "residuals"))
In [59]:
prediction_real(regression_taxi.taxi_call_num.values, models[80].resid.values, (0, 168), "first january week", 'residuals',
(15, 5), ("real values", "residuals"))
not accounted structure stay in residual plot¶
In [61]:
plt.figure(figsize(80,10))
sm.tsa.seasonal_decompose(models[80].resid).plot()
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(models[80].resid)[1])
In [92]:
#sm.tsa.seasonal_decompose(models[80].resid).__dict__.keys()
#plt.plot(sm.tsa.seasonal_decompose(models[80].resid).__dict__["resid"][:96])
In [62]:
regression_taxi.iloc[:, :321].head()
Out[62]:
In [66]:
full_regres_taxi = regression_taxi.iloc[:, :321]
full_regres_taxi.shape
Out[66]:
its possible to introduce additional feature to consider snow storm in january 26-27. But this will not help us in making future prediction.¶
In [67]:
#full_regres_taxi = full_regres_taxi.drop(['jan_anom'], axis=1)
full_regres_taxi['jan_anom'] = [full_regres_taxi.loc[date].taxi_call_num if date.month == 1 and date.day in [27] and date.year == 2015 else 0 for date in full_regres_taxi.index]
full_regres_taxi.shape
Out[67]:
In [69]:
m1 = smf.ols('taxi_call_num ~ ' + ' + '.join(list(full_regres_taxi.columns)[1:]) + ' + jan_anom', data=full_regres_taxi)
fin_model_anomal = m1.fit(cov_type='HC1')
In [70]:
#result = prediction(fin_model_anomal, components=81, feature_num=4, dataframe=regression_taxi)
coefs = np.array([fin_model_anomal.params[coef] for coef in fin_model_anomal.params.index])
fourier_components = np.hstack((np.array([1]*full_regres_taxi.shape[0]).reshape(full_regres_taxi.shape[0],1),
full_regres_taxi.iloc[:, 1:4*80+2].values))
result = coefs.dot(fourier_components.T)
In [71]:
prediction_real(regression_taxi.taxi_call_num.values, result, (744, 1416), "february", 'taxi calls',
(15, 5), ("real values", "prediction"))
Now we are going to select appropriative ARIMA model for residual signal - data not encountered by linear regression model¶
In [74]:
def plot_autocorr(resid_data):
plt.figure(figsize(15,10))
ax = plt.subplot(211)
sm.graphics.tsa.plot_acf(resid_data.squeeze(), lags=168, ax=ax)
pylab.show()
ax = plt.subplot(212)
sm.graphics.tsa.plot_pacf(resid_data.squeeze(), lags=168, ax=ax)
pylab.show()
In [75]:
plot_autocorr(fin_model_anomal.resid.values)
single ordinary differentiation of the time series¶
In [77]:
full_regres_taxi["residuals"] = fin_model_anomal.resid
full_regres_taxi["residual_diff1"] = full_regres_taxi["residuals"] - full_regres_taxi["residuals"].shift(1)
In [78]:
plot_autocorr(full_regres_taxi["residual_diff1"][1:].values)
select search parameters from autocorrelation and partial autocorrelation plots (autocorrelation - q, Q)¶
In [41]:
ps = range(0, 12)
qs = range(0, 6)
Ps = range(0, 2)
Qs = range(0, 2)
In [42]:
parameters = product(ps, qs, Ps, Qs)
parameters_list = list(parameters)
len(parameters_list)
Out[42]:
In [52]:
%%time
residuals_results = []
residuals_best_aic = float("inf")
warnings.filterwarnings('ignore')
resid_wrong_param = 0
for counter, param in enumerate(parameters_list):
if counter % 10 == 0:
print(counter, end=' ')
#some of the parameters set are incompatible
try:
model=sm.tsa.statespace.SARIMAX(full_regres_taxi["residuals"], order=(param[0], 1, param[1]),
seasonal_order=(param[2], 0, param[3], 24)).fit(disp=-1)
#print failure parameters and continue
except LinAlgError:
print("singular matrix: ", param)
continue
except ValueError:
resid_wrong_param += 1
continue
aic = model.aic
#save best model, aic, parameters
if aic < residuals_best_aic:
residuals_best_model = model
residuals_best_aic = aic
residuals_best_param = param
residuals_results.append([param, model.aic])
warnings.filterwarnings('default')
In [53]:
print(sorted(residuals_results, key=lambda tup: tup[1], reverse=False)[:5])
In [54]:
plt.figure(figsize(15,12))
ax = plt.subplot(211)
sm.graphics.tsa.plot_acf(residuals_best_model.resid.values.squeeze(), lags=168, ax=ax)
pylab.show()
ax = plt.subplot(212)
sm.graphics.tsa.plot_pacf(residuals_best_model.resid.values.squeeze(), lags=168, ax=ax)
pylab.show()
residual signal from best selected ARIMA model looks completely decorrelated¶
In [79]:
best_arima = sm.tsa.statespace.SARIMAX(full_regres_taxi["residuals"], order=(11, 1, 1), seasonal_order=(1, 0, 1, 24)).fit(disp=-1)
In [55]:
plt.figure(figsize(15,8))
plt.subplot(211)
residuals_best_model.resid[1:].plot()
plt.ylabel(u'Residuals')
ax = plt.subplot(212)
sm.graphics.tsa.plot_acf(residuals_best_model.resid[1:].values.squeeze(), lags=168, ax=ax)
print("Критерий Стьюдента: p=%f" % stats.ttest_1samp(residuals_best_model.resid[1:], 0)[1])
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(residuals_best_model.resid[1:])[1])
In [80]:
prediction_real(full_regres_taxi["residuals"], best_arima.fittedvalues, (744, 1416), "february", 'taxi calls',
(15, 5), ("real residuals", "arima fitted values"))
In [81]:
prediction_real(full_regres_taxi["residuals"], best_arima.fittedvalues, (3624, 4344), "june", 'taxi calls',
(15, 5), ("real residuals", "arima fitted values"))
ARIMA fitted values looks good¶
In [82]:
best_arima.fittedvalues.shape
Out[82]:
add ARIMA with parameters (p, q, P, Q) = (11, 1, 1, 1) fitted values to values predicted by linear regression model and get final result¶
In [83]:
final_res = result + best_arima.fittedvalues
plot final fit for february and june mounth¶
In [85]:
prediction_real(regression_taxi.taxi_call_num.values, final_res.values, (3624, 4344), "june", 'taxi calls',
(15, 5), ("real taxi calls", "final prediction"))
In [86]:
prediction_real(regression_taxi.taxi_call_num.values, final_res.values, (744, 1416), "february", 'taxi calls',
(15, 5), ("real taxi calls", "final prediction"))
In [ ]:
Try to fit data using only ARIMA model without linear regression on fourier features¶
look at ACF и PACF of time serie:¶
In [6]:
diff_data=pd.DataFrame(taxi_data.taxi_call_num, columns=["taxi_data.taxi_call_num"])
In [88]:
plot_autocorr(taxi_data.taxi_call_num.values)
performe ordinary differentiation¶
In [7]:
diff_data["diff_168"] = taxi_data.taxi_call_num - taxi_data.taxi_call_num.shift(168)
In [90]:
plot_autocorr(diff_data["diff_168"][168:].values)
performe seasonal differentiation with weekly period (168 hours)¶
In [8]:
diff_data["diff_168_1"] = diff_data["diff_168"][168:] - diff_data["diff_168"][168:].shift(1)
In [92]:
plot_autocorr(diff_data["diff_168_1"][169:].values)
In [26]:
ps = range(0, 12)
qs = range(0, 6)
Ps = range(0, 2)
Qs = range(0, 2)
In [27]:
parameters = product(ps, qs, Ps, Qs)
parameters_list = list(parameters)
len(parameters_list)
Out[27]:
In [28]:
%%time
results = []
best_aic = float("inf")
warnings.filterwarnings('ignore')
wrong_parameters_counter = 0
for counter, param in enumerate(parameters_list):
if counter % 10 == 0:
print(counter, end=' ')
#some of the parameters set are incompatible
try:
model=sm.tsa.statespace.SARIMAX(taxi_data.taxi_call_num, order=(param[0], 1, param[1]),
seasonal_order=(param[2], 1, param[3], 24)).fit(disp=-1)
#print failure parameters and continue
except LinAlgError:
print("singular matrix: ", param)
continue
except ValueError:
wrong_parameters_counter += 1
continue
aic = model.aic
#save best model, aic, parameters
if aic < best_aic:
best_model = model
best_aic = aic
best_param = param
results.append([param, model.aic])
warnings.filterwarnings('default')
In [29]:
print(wrong_parameters_counter)
In [30]:
print(sorted(results, key=lambda tup: tup[1], reverse=False)[:5])
look at residuals of best selected SARIMA model¶
In [31]:
plt.figure(figsize(15,12))
ax = plt.subplot(211)
sm.graphics.tsa.plot_acf(best_model.resid.values.squeeze(), lags=168, ax=ax)
pylab.show()
ax = plt.subplot(212)
sm.graphics.tsa.plot_pacf(best_model.resid.values.squeeze(), lags=168, ax=ax)
pylab.show()
In [32]:
plt.figure(figsize(15,8))
plt.subplot(211)
best_model.resid[1:].plot()
plt.ylabel(u'Residuals')
ax = plt.subplot(212)
sm.graphics.tsa.plot_acf(best_model.resid[1:].values.squeeze(), lags=168, ax=ax)
print("Критерий Стьюдента: p=%f" % stats.ttest_1samp(best_model.resid[1:], 0)[1])
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(best_model.resid[1:])[1])
Согласно критерию Стьюдента - остатки не смещены
In [9]:
best_arima_1 = sm.tsa.statespace.SARIMAX(taxi_data.taxi_call_num, order=(9, 1, 1),
seasonal_order=(1, 1, 1, 24)).fit(disp=-1)
plot final fit for february and june mounth¶
In [15]:
prediction_real(regression_taxi.taxi_call_num.values, best_arima_1.fittedvalues.values, (3624, 4344), "june", 'taxi calls',
(15, 5), ("real taxi calls", "final SARIMA fitted"))
In [16]:
prediction_real(regression_taxi.taxi_call_num.values, best_arima_1.fittedvalues.values, (744, 1416), "february", 'taxi calls',
(15, 5), ("real taxi calls", "final SARIMA fitted"))
Using only SARIMA (parameters (p, q, P, Q) = (9, 1, 1, 1)) prediction model with one ordinary and one seasonal (168 hour) differentiation we got visualy good fitted values for time serie. SARIMA is very time and resource consuming especially with seasonal differentiation. It give good fitted results but is pure for prediction of such taxi trips time serie.¶
In [ ]: